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Abstract. Using a so-called hemispherical model we derive a general transport equation for cosmic ray and 
thermal particles scattered in pitch angle by magnetic inhomogeneities in a moving collisionless plasma. The weak 
scattering through 90 degrees results in isotropic particle distributions in each hemisphere. The consideration 
is not limited by small anisotropies and by the condition that particle velocities are higher than characteristic 
flow velocity differences. For high velocities and small anisotropies the standard cosmic ray transport equation is 
recovered. We apply the equations derived for investigation of injection and acceleration of particles by collisionless 
shocks. 
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1. Introduction 

Progress in cosmic ray astrophysics has been made with 
the introduction of the transport equation for cosmic rays 
(KrymskyQMl Parker HMSl Dolginov & Toptygin 
Gleeson & Axford I1969p . It was successfully applied to 
many astrophysical problems: cosmic ray modulation in 
the solar wind, acceleration of particles at shock fronts, 
cosmic ray propagation in the Galaxy etc. This equation 
implies that in the strong scattering approximation the 
particle distribution is almost isotropic, which in partic- 
ular implies that velocities of particles are higher than 
the characteristic flow velocity differences in the plasma. 
However, in several circumstances these conditions are vi- 
olated and more general kinetic equations should be used 
which also deal with the evolution of the pitch angle dis- 
tribution. Since the solution of such equations is not a 
simple task, we would like find a description that comes 
close to that implied by the standard transport equation. 

Energetic particles are scattered in pitch angle relative 
to the mean magnetic field direction by random magnetic 
inhomogeneities. It is well known that in the case of par- 
ticle transport along the mean magnetic field a so-called 
problem of scattering through the pitch-angle of 90 de- 
grees exist. According to quasilinear theory the resonant 
scattering of particles is weak at small pitch-angle cosines 
(Jokipii I1966|) . The main reason for this phenomenon is 
the small energy density of short waves propagating along 
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the magnetic field and scattering the particles with pitch 
angles close to 90 degrees. These waves may be easily 
damped by thermal ions. Even without such a damping, 
the scattering in the vicinity of 90 degrees is depressed 
in the case of a power-law spectrum of waves that is pro- 
portional to fc~ 2 or steeper. Here k is the wavenumber. 
This problem can be avoided if one takes into account 
the scattering by oblique waves or nonlinear interactions, 
for example magnetic mirroring (see e.g. V51k ll975l 119751 
Achterberg |198ip . One can expect that the scattering ef- 
ficiency is small near the pitch angle of 90 degrees. 

Observations of suprathermal pick-up particles in the 
solar wind suggest that this is the case (Fisk et al. I1997p . 
If the scattering through 90 degrees is weak, the parti- 
cle distribution is almost isotropic in each hemisphere of 
the velocity directions parallel or antiparallel to the mean 
field. So, in this case the angular dependence of particle 
distribution is reduced to the two number densities of for- 
ward and backward moving particles. The corresponding 
equations were derived by Isenberg (|1997|) for the case of 
pick-up ions in the solar wind. 

We use this approach and derive general transport 
equations for arbitrary nonrelativistic flow velocities. We 
use the particle distributions in the frame moving with 
the medium. This permits us to consider the case of arbi- 
trary particle velocities. Since the number density of par- 
ticles in both hemispheres may substantially differ, large 
anisotropies can also be taken into account. The equations 
derived have a broad range of applicability: propagation 
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of solar energetic particles and pick-up ions in the solar 
wind, injection into diffusive shock acceleration and other 
processes. 

As an example of the application of the equations de- 
rived, we consider the problem of diffusive shock accelera- 
tion. Since our equations describe propagation of thermal 
as well as energetic particles in the vicinity of the collision- 
less shock front, they can be used for the determination of 
the shock velocity profile and the spectra of thermal and 
nonthermal particles without any assumptions regarding 
the injection efficiency. One only has to prescribe the scat- 
tering law. In this sense the approach is similar to the one 
used by Ellison and Eichler (|1984p . They applied a Monte- 
Carlo technique for the solution of the kinetic equation. 

The basis of such an approach is the idea that the in- 
teraction of upstream and downstream plasmas may result 
in instabilities similar to firehose instability (Parker 1961, 
Quest I1988|) . The magnetic inhomogeneities produced by 
such instabilities may play the role of scattering centers 
and provide the shock dissipation and heating. These ideas 
are also the basis of the so-called hybrid simulations of col- 
lisionless shocks (Leroy & Winske[1983 Quest [19861) . Th e 
ions move in self-consistent electromagnetic fields and the 
electrons are considered as a charged neutralizing fluid in 
these simulations. 

The equations are derived in the next two Sections. 
Application to the case of collisionless shock is considered 
in Sect. 4 and 5. Sect. 6 contains a discussion of results 
obtained and conclusions. 



2. Basic equations 

We start with the kinetic equation for the momentum dis- 
tribution /(r,p) of charged particles 



^+vV/ + F i |^ = 5i/ 
at op 



(1) 



Here v and p are the particle velocity and momentum 
respectively, Fx, = q~E + (q/c)[v x B] is the Lorence force 
and the operator St describes the pitch angle scattering of 
particles by random magnetic inhomogeneities in the fluid 
frame. We shall further assume that the mean electric field 
E can be described as 



1, 



E = Eub - -[u x B] 



(2) 



where b is the unit vector along the direction of the mean 
magnetic field B, u is the mass velocity and E\\ is the 
component of the mean electric field which is parallel to 
the mean magnetic filed. 

It is convenient to perform a change of variable p to 
p' that formally coincides with the nonrelativistic Lorence 
transform 



(3) 



Averaging the particle gyrophase angle and neglecting 
terms of the order u/c << 1 we find (cf. Skilling [19711 
KulsrudQSM]): 



% + («Vb + u)V/ = Stf- 
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(4) 



Here // is the cosine of the pitch angle of the particle and 
F\\ is the sum of the parallel components of the electric 
and inertia force: 



p' ( du 

^ll=^||-bM^ + (uV)u 



The scattering operator can be written as 
- d 1 - // 2 , df 



(5) 



(6) 



Here v is the scattering frequency. 

Let us assume that the scattering frequency is large for 
pitch angle cosines > fiQ. The momentum distribution 
/ is then isotropic in these regions of the angle space. We 
introduce distributions of particles in the backward and 
forward hemispheres N±: /(p', //) = N-(p') for (j! < — /i 
and f(p',fJ,') — N + (p') for p! > fi . Integration of Eq. (4) 
from -1 to — fiQ and from fj,Q to 1 gives the equations for 
these distributions: 



dN± 



vdf_ 
2 dpi 1 



u± — b 

2 



=±A»o 



1 p' \ dN± 
±-F|i - — Vu — — 

2 11 3 J dp' 



(7) 



It was assumed that /io < < 1 . Now we calculate the right- 
hand side of these equations. For this we must find the 
solution of Eq. (4) in the region |//| < ^o- Since fj,o << L 
we can leave only the terms with derivatives on // in this 
equation with the result: 



p 1 J dfi' d/i' 2 d/i' 



(8) 



This equation should be solved with the boundary con- 
ditions /(p',±/io) = N±(p'). Substitution of the solution 
into the Eq. (7) then gives 

dN± 



at 



u±-b 



4 F n 



-Vu 



dN± 
dp' 



(9) 



Here the frequencies v± describe the rate of particle ex- 
change between forward and backward hemispheres: 



v±(p')=±- 



exp 



±A 



\v'p' 



Vb 



(10) 
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and the mean free path A of the particle in the uniform 
medium is given by the expression: 



A 



(11) 



A very similar equation was derived by Isenberg (1997) 
for pick-up ions in the solar wind. We found more accurate 
expressions for the frequencies v± and these coincide with 
the result of Kota (|2000|) obtained for the case Fii = 0. 

We neglected the derivatives in time and space in Eq. 
(8). This assumption is valid if the characteristic time r 
and the characteristic scale I of the problem considered are 
large enough: r, l/v' >> ^/v ~ u$\/v' . Since fi « 1, 
our derivation is valid even if the scale I is comparable 
with the mean free path A. The equation (9) is exact in 
the mathematical limit /iq — > 0, A — » const. 

The Eqs. (9) can be rewritten in the conservative form: 



u± — b 

2 



dN± 
1 9 , 2 / 1 



Vu N ± 



v±N±. (12) 



The equations derived have a simple physical meaning. 
The second term in the left hand side describes the trans- 
port of the particles by the medium moving with velocity 
u and the proper motion of particles along the magnetic 
field with the speed ±u'/2 (to be compared with "coher- 
ent" propagation of solar energetic particles considered by 
Earl (|1974p ). The third term describes the energy losses 
or gains in the inhomogeneous flow and electric field. The 
right hand side corresponds to the exchange of particles 
between forward and backward hemispheres. 

It is clear that the system of equations (12) is hyper- 
bolic. It can be reduced to uncoupled equations for N + 
and A_ : 



dN± 



p'dN±^ 



J__d_ ,-. 

p/2 Qp/P 



3 dp' 
P' 



d_ 

di 



Vu 



P „ \ dN ± 
-Vu 



dp' 



Another interesting case is when the energy changes 
and transport by the medium are negligible. Then Eq. (13) 
is reduced to the so-called telegraph equation (cf. Fisk & 
AxfordCEHMl): 



dN± 

~dT 



dt 



-(Vb) 



l 

2iZ 



dN± 



± -(bV)JV± 



.(15) 



This equation describes propagation of cosmic ray par- 
ticles with finite velocity and can be reduced to the diffu- 
sion equation only in the case of a slow time evolution. 

3. Self-consistent closure 

The equations derived in the previous section contain pre- 
scribed values for the flow velocities u, the parallel force 
.Fjl and the magnetic field B. Such a description is a good 
approximation for the investigation of propagation of test 
particles. However, in many interesting cases these quan- 
tities should be determined self-consistently, involving the 
distributions N + and i\T_. 

We treat electrons as a massless fluid. In this case the 
Euler equation for electrons is reduced to an expression 
for the electric field: 



V P 1 

E + SE = [u e x (B + SB)], 

qn e c 



(16) 



where u e and P e are the velocity and the pressure of the 
electron fluid respectively, n e is the electron number den- 
sity, SE and 5H are the random components of the electric 
and the magnetic fields respectively. The velocity u e can 
be found from Maxwell's equation 



[V x (B + SB) 



c 



d 3 p(v-u e )(/ + ( 5/), 



(17) 



where Sf is the random component of the ion distribution, 
and we have neglected the displacement current by assum- 
ing a slow time evolution. It was also assumed that the 
plasma is pure hydrogen and quasineutral. Substitution of 
u e from this equation into Eq. (16) and averaging gives: 



E=--[uxB]-- 

c qn % 



[B x [V x B]] + VP e 



.(13) 



d 3 p P Stf 



(18) 



In the case of high scattering frequencies the parti- 
cle distributions in the different hemispheres are almost 
equal to each other, N + « N- w N. Assuming also a slow 
time evolution and particle velocities much larger than 
the medium velocity we come to the standard transport 
equation for cosmic rays (Krymsky I1964[ Parker [1965, 
Dolginov & ToptyginElSl Gleeson & Axford [MSI : 

dN n' dN 

l^+uViV-I^Vu^ (Vb)£>,|(bV)JV. (14) 

Here the diffusion coefficient along the magnetic field 
D\\ = u'A/4. 



Here rij is the mean number density of ions and we 
have used the formal definition of the scattering opera- 
tor Stf = — 2 < jp[(v — u) x SB]Sf >, where the angle 
brackets mean the average over the random fluctuations of 
SB and Sf of the magnetic field and the momentum distri- 
bution. We also neglected the magnetic tension of the ran- 
dom component. Generally speaking, the scattering oper- 
ator does not conserve momentum. Since the magnetic 
inhomogeneities are frozen into the electron fluid, this ad- 
ditional force acts on the electron fluid. As a result, an 
additional electric field due to the charge separation ap- 
pears. It is described by the last term in the parentheses 
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of Eq. (18). Such a field is indeed observed in hybrid sim- 
ulations of collisionless shocks (Quest [T"988|) . 

We multiply Eq.(l) containing the electric field (18) 
by pd 3 p and integrate over momentum space. This gives 
the Euler equation of motion: 



| + (uV)u)=-i[Bx|VxB] 



Here 

Pi = l f d 3 p'v'p'(N + + N-) 



V(P e + Pi)-(19) 



(20) 



is the pressure of ions and p is the mean density. The 
parallel electric field can be found from Eq.(18). 

Since the Eq. (1) with the electric field (18) conserves 
momentum, the same is true for Eqs. (9). Let us multiply 
them by p' ', integrate over p ,2 dp' and subtract the first 
from the second. This gives after some algebra 



1 

Tii 

du 

~dt ' 



d 6 p'p'(v + +^)(N+ 



bViV 



(JV++JV-) 



(21) 



The last term term in the right hand side of this expres- 
sion can be neglected, since the force F\\ is essential only 
for nonrelativistic particles which determine the plasma 
density. 

We shall assume an adiabatic evolution of the electron 
pressure 

dP, 



it + uVP e + -P e Vu = 0, 
at 3 



(22) 



and frozen-in magnetic field. The last term in Eq.(18) is 
essential for momentum conservation but can be neglected 
in Faraday's induction equation: 

<9B 



at 



= [V x [u x B]]. 



(23) 



We have thus obtained the closed system of Eqs. (12), 
(19), (22), (23) with expressions (10) for the frequencies 
v± and Eq. (21) for the parallel force P||. In the next sec- 
tions we apply these equations to the combined problem 
of injection and nonlinear acceleration at a plane shock. 

4. Application for acceleration at the plane shock 

Let us consider the acceleration at the parallel one- 
dimensional shock. We can write the steady-state version 
of Eq. (9) for nonrelativistic particles in the shock frame 
as 



tA dN± f F\ v'du\ dN ± 
U 2 dz + { 2m~ 3 dz dv' 



(24) 



Here m is the mass of particles and the parallel force Fn = 
qE\\ — mudu/dz. 



0.5 




Fig. 1. Characteristics of the hyperbolic system of equa- 
tions (24) for the case E\\ = 0. Characteristics for forward 
moving particles and backward moving particles are shown 
by solid and dotted lines respectively. 

The characteristics of this system of equations for the 
case P|| = are given by 



u 3 ~^ 3 



V? 



u 3 3 



const (25) 



They are shown in Fig.l. An upstream particle begins 
its motion with low velocity from the upper left part of 
this u — v 1 plane and goes down the characteristics. The 
forward moving particles gain energy, backward moving 
ones lose energy. At any moment the particle may change 
hemisphere and will continue the motion along another set 
of characteristics. If the compression ratio of the shock 
is high enough, the forward moving particle can change 
hemisphere in the vicinity of the line u = i>'/2 (shown by 
the dash- dotted line) and return upstream along the back- 
ward characteristic with energy gain. It can again change 
the hemisphere upstream and again move in a downstream 
direction etc. This may be repeated many times and the 
accelerated particle goes to the right beyond the part of 
u — v' plane shown in Fig.l. 

There is an another possibility for the initial acceler- 
ation. If the initial speed of the backward moving parti- 
cle upstream is high enough, say 0.7wi, it can reach the 
line u = v' jl directly and return upstream with an en- 
ergy gain. However, the initial speed is rather high. The 
number of such particles is small for high Mach number 
shocks. 

The characteristics of Eq. (24) depend on Fn . The char- 
acteristics for the case P|| = are shown in Fig. 2. This case 
is close to the real situation because the two first terms 
in Eq. (21) almost cancel each other (see also the numeric 
modeling below). In other words, the electric force qEu al- 
most compensate the inertia force —mudu/dz in the sec- 
ond term on the left-hand side of Eq. (24). This is not a 
simple coincidence. The electric force is the consequence of 
the momentum transfer and is approximately —dP/dz/p 



V.N. Zirakashvili: Hemispherical transport equation 



5 
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u/u. u = u 1 



u=v'/V 




Fig. 2. Characteristics of the hyperbolic system of equa- 
tions (24) for the case F\\ = 0. Characteristics for forward 
moving particles and backward moving particles are shown 
by solid and dotted lines respectively. 



that is exactly udu/dz. The characteristics are determined 
by expressions 



v' 4 ± —uv' 3 — const. 



(26) 



We find the analytical solution of Eq. (24) for the case 
when the velocity profile is discontinuous: u{z) = iti, z < 
and u(z) — u 2 , z > 0. In the regions of constant velocity u 
the system (24) can be rewritten as (to be compared with 
Gombosi et al. [1993]) 



8N+ d 



= — (J* - 4u z 



a on+ 



dz dz Av' dz 

N_=N+ + (2u + v') ~ 



(27) 
(28) 



XdN+ 

' v 1 dz ' 

It was also assumed that En = 0. The solution of these 
equations upstream (z < 0) is 

N ± = Noo + (jv 1± - JVoo) exp 4 l lV ' Z / X 2 . (29) 

t/ z — 4m| 

Here N^, is the distribution at z — — oo, Ni + and Ni- 
are the distributions N + and iV_ just upstream the shock. 
According to Eq. (28) they are related as 

(«' + 2ui)(JVi+ - TVoo) - («' - 2ui)(JVi_ - N^, v' > 2u x 

N 1+ = Ni- =N OD , v'< 2u x (30) 

The solutions downstream are given by 



N ± =N 2+ [- + — ) + N 2 -(-- — 
" \2 4u 2 ) \2 Au 2 



Here N 2+ and N 2 - are the distributions just downstream. 
They are related as 



N 2 + = 2V 2 -, v' > 2u 2 . 



(32) 



This means that distributions N + and AL. do not depend 
on z for v' > 2u 2 . 

We find the relation between the upstream and down- 
stream distributions N\± and N 2 ±. It depends on the so- 
lution of Eq. (24) in the transition region. We consider 
the case Fn =0. Then the distributions upstream and 
downstream are related by characteristics (26): 

N 1± = F ± (v' 4 ± ^uW 3 ), N 2± = F ± (v' 4 ± ^u 2 v' 3 ). (33) 

Here F± are two functions to be determined. Thus we 
have the six equations (30), (32), (33) for the six unknown 
functions Nx±, N 2 ± and F±. 

Let us consider the case N^ = 6(v' — vq)/(2Vq). Since 
the change of velocity in the transition region is governed 
by the characteristics (26), the solution is the following 
sum of 5-functions: 

2 00 

N 1± =A 8{v l -v,)/{v f+Y J Y.^± 5{ - v '-<)l { -O 2 ^ M ) 

3 = 1 i=l 



2 00 
j=l i=0 



(35) 



Here A 3 i± , B 3 i± , v\ and w 3 are sequences to be determined. 
The index j may have the values 1 and 2. These values 
correspond to the two possibilities of injection of forward 
and backward moving particles. They were described after 
Eq. (25). The different values of the index i correspond 
to consequent states of the one particle moving between 
upstream and downstream regions of the shock. 

If the upstream sequences Aj± are known for v' > 2u\ 
the downstream ones can be found from characteristics 
(26). Using the properties of 5-functions we find that they 
are given by 



B{ ± = V { + 2Ul A %+ , w{ > 2u 2 



w\ + 2u 2 



(36) 



and the downstream velocity w? is the solution of equation 



(«l) 4 + ^iW) 3 = K) 4 + ^2K) 3 - 



(37) 



The next Aj +1± are 



wj — 2u 2 



B 



id - 2u 2 



v i+i - 2u i ' v i+i + 2u i 

where the velocity v? , 1 is the solution of the equation 



(38) 



. , . 4u 2 v'z/X 
^~ 2 ){N 2 --N 2+ )^-i-^. 



(31) 
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Eqs. (36)-(39) can be used recurrently to calculate the 
sequences A\ ± , B{ ± , v\ and w\ . 

Let us introduce the critical velocities v c ±. They are 
the solutions of equations: 



vl±(3v c ± ± 8wi) = (48 ± 64)w|. 



(40) 



The initial velocity v and A a = 1/2 are given. There 
are four cases depending on the initial velocity v and 
velocities v c ±. 

1) vo < v c -. Forward and backward moving particles 
reach the downstream region with velocities smaller than 
2u2 and are further advected. The injection in acceleration 
does not occur. The downstream coefficients Bq± are given 
by 



B o+ 



A 



v + 2«i 
-2u 2 

2m 



Bl 







Wn 



Wn 



2u 2 



B l _ = 



(41) 



10 r 



0.1 



0.01 



N(v')V 



Fig. 3. The analytic solution for u\ju 2 = 4 and vo = 
0.25ui for the case F\\ = 0. The downstream and upstream 
distributions are shown by solid and dashed lines respec- 
tively The height of the vertical lines corresponds to the 
coefficients in Eqs. (34) and (35) 



where the velocities and Wq are the solutions of equa- 
tions: 



and the downstream velocity u>J is the solution of the 
equation 



(v ) 4 + ^u 1 (v f = (w 1 ) i + -u 2 (w 1 Q f 



(«o) 4 +3«i(«o) 3 



(46) 



M 4 3«iM 3 



2\4 



K) 



3"2(^o) 3 - 



(42) 



All other coefficients A\ ± = 0, Bf ± = for i > 0. 

2)v c - < v o < v c+ . The most important case for injec- 
tion. The forward moving particle reaches the downstream 
region with a velocity less than 2u 2 and again is not in- 
jected into acceleration. The coefficients B}± and A\ ± are 
the same as in the previous case. 

The backward moving particle goes along characteris- 
tics and returns upstream (see Fig. 2). Using properties of 
(5-functions in Eq. (33) and conditions (30) we found 



A\_ = A 



vq — 2u\ 



2ui 



A 2 



At 



2ui 



2ui 



Here the velocity v\ is the solution of the equation 



2\3 



(«0) 4 -3«l(«0) 3 = («l 2 ) 4 -g«lK) 



(43) 



(44) 



The coefficients A\ ± can be found from equations Eqs. 
(38) and (39). Now the coefficients B}± at i > and Aj ± 
at i > 1 can be found consecutively using Eqs. (36)- (39). 

4) v > 2u\. The coefficients A J i± and Bf ± at i > can 
be found consecutively using Eqs. (36)- (39). 

Using Eqs. (36) and (38) one can relate Aj +1+ and 



A\ + - 



if, 



A 1 - A 1 



2u 2 v 1 ; + 2ui 



vj +1 + 2u\ w\ + 2u 2 



(47) 



It is easy to show that for large wj and vf the change 
of velocity 



-(til -u 2 ). 



+i - " , ~ - "2i- (48) 

Using Eq. (47) we find that for high velocities 

(49) 



It is clear that vf is close to 8ui/3 (see Fig. 2). The other 
coefficients A 2 ± and Bf ± can be found consecutively using 
Eqs. (36)-(39). 

3) v c + < vo < 2u\. Forward and backward moving par- 
ticles reach the downstream region with velocities greater 
than 2u 2 and are injected into acceleration. The coeffi- 
cients Af ± and Bf± are the same as in the previous case. 
The coefficients Bq± are given by 



Di vq + 2 Ui 



Wn 



2u 2 



(45) 



that is exactly the spectrum of particles produced by dif- 
fusive shock acceleration. 

The solution for u\/u 2 — 4 and v — 0.25ui is shown 
in Fig. 3. The critical velocities are v c _ w 0.2ti! and 
v c+ w 0.35ui in this case. The far upstream distribution is 
shown by the dashed line on the left. The forward moving 
particles of this distribution reach the downstream region 
with slightly higher velocity (solid line on the left). The 
backward moving particles are directly accelerated along 
the characteristic and return upstream with a velocity 
close to 8«i/3 (first dashed line on the right). Particles are 
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I ' ' ' ' 1 ' ' ' ' 1 ' ' ' ' V - I - - , - - - , - -r - - , - - 

-1 -0.5 0.5 1 

Fig. 4. The results of simulation of the parallel shock with 
Mach number 3.77, T e /T; = and energy independent 
free path A = 0.2. The forward and backward particle 
distributions downstream N + and iV_ are shown on the 
top panel by solid and dashed lines respectively. Plasma 
velocity u, pressure of ions P and electric force F = qEu 
are shown on the bottom panel by solid , dashed and dotted 
lines respectively. The total compression ratio r = 4.05 is 
obtained. 
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Fig. 5. The results of simulation of the parallel shock with 
Mach number 7.75, T e /Ti = and energy independent 
free path A = 0.2. The forward and backward particle 
distributions downstream N + and iV_ are shown on the 
top panel by solid and dashed lines respectively. Plasma 
velocity u, pressure of ions P and electric force F = qEu 
are shown on the bottom panel by solid , dashed and dotted 
lines respectively. The total compression ratio r = 6.67 is 
obtained. 



accelerated further, moving between downstream and up- 
stream regions of the shock (the other lines on the right). 
The similar (S-function spectrum was obtained by Bogdan 
and Webb {1987) for cosmic ray acceleration in the so- 
called two-streaming approximation of Fisk and Axford 

The spectrum will be smoother if one takes into ac- 
count the scattering of particles in the transition region. 

5. Numeric modeling of collisionless shocks 

We use the system of equations (12), (19), (22), (23) for the 
modeling of steady-state one-dimensional nonrelativistic 



shocks. We neglect the dynamical effects of the magnetic 
field, which means that the magnetic field produces only 
kinematic effects that are essential for oblique shocks. The 
upstream plasma state is then determined only by the 
sonic shock Mach number M and by the ratio of electron 
and ion temperatures T e /Ti. 

A so-called Alfven heating is not included in our 
model. It is well known that the energetic particles up- 
stream of the shock can effectively generate Alfven waves 
due to streaming instability (Wentzcl 1969] Bell I1978p . 
The damping of these waves results in the heating of the 
shock precursor (Volk & McKenzie 119821 McKenzie & 
Volk ll982"[) . This heating is essential, since it decreases the 
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effective Mach number and compression ratio of the shock. 
As shown in the Appendix, this effect can be formally in- 
cluded using the relation between sonic and Alfven Mach 
numbers of the shock: 



M 



MiM„ 



M a + §Mf 



(50) 



The results obtained below for shocks with sonic Mach 
number M without Alfven heating can be formally used 
for shocks with Alfven heating and with Alfven and sonic 
Mach numbers M a and M s respectively. 

The details of the numeric method are given in the 
Appendix. The plasma flow with unity velocity and unity 
density enters the simulation box from its left boundary. 
The pressure at the right boundary was adjusted accord- 
ing to the motion of the shock in order to reach the steady 
state. We used the free escape boundary condition for the 
particle distribution at the left boundary. 

We prescribe the dependence of the free path of parti- 
cles on the velocity v' and space coordinate z. We used the 
free path A that is independent on z. This is not a strong 
limitation. Indeed, in the rather general case of separable 
dependence of A on v' and z, the results obtained depend 
only on the integral J dz/X(z). Thus our results can be 
used also for this more general case. 

The energy dependence of the mean free path A is de- 
termined by the spectrum of magnetic inhomogeneities 
and by the mechanism of scattering in the vicinity of a 
pitch angle of 90 degrees. Since the low frequency waves 
propagating along the magnetic field are subject to non- 
linear wave steepening, one can expect that the magnetic 
spectra are close to k~ 2 and the corresponding free path 
does not depend on energy. Such spectra (or slightly flatter 
ones) were indeed observed in the hybrid plasma simula- 
tions of parallel collisionlcss shocks (Giacalone ct al. 1993, 
Giacalone 2004) and in the vicinity of the Earth bow shock 
in the solar wind. In the results presented here we used 
the energy independent free path. 

We limit ourselves to the case of cold electrons T e /Ti = 
0, as the dissipation of the Alfven waves in the precursor 
may result in the preferable heating of ions. The solar 
corona is an example of this. 

We start with the case of parallel shocks with Mach 
numbers 3.87 and 7.75. According to Eq. (50) the corre- 
sponding Alfven Mach numbers are 7.5 and 30 that can be 
applied to the interplanetary and supernova shocks respec- 
tively. The downstream particle distributions are shown in 
the top panels of Fig. 4 and 5. The velocity and gas pres- 
sure profiles are shown in the bottom panels. 

The simulated shocks demonstrate clear sub-shocks 
with a width of about one half of the free path A and 
a smooth precursor produced by the particles in the high 
energy tail that is seen in the top panels of Fig. 4 and 5. 
It was found that the main part of the downstream distri- 
bution in Fig. 4 and 5 is the adiabatically compressed far 
upstream distribution which was taken to be Maxwellian. 
The total electric potential corresponding to the electric 
force qE\\ is about 0.2-0.3 mu\ according to Fig. 4 and 5. 



The total compression ratio of the shock can be found 
from the following simple estimate. 

In the precursor region the thermal particles are heated 
adiabatically. When their velocities become comparable to 
the plasma velocity, the backward moving particles are ac- 
celerated more efficiently (see Fig. 2) and injected for fur- 
ther acceleration. The compression ratio should be large 
enough for this to be the case. For estimation we assume 
that distribution of upstream particles is proportional to 
5(v' — vt), where vt = ^/3/5ui/M is the thermal velocity 
of the plasma. The maximum shock downstream velocity 
ui that is necessary for injection of such particles accord- 
ing to Fig. 2 and characteristics (26) is 



u 2 



1 U Ul 

-VT 8 — 

2 \ v T 



1/4 



(51) 



Thus the thermal velocity vt should be larger than the 
critical velocity u c _ (see Sect. 3). The second term in 
parentheses can be neglected. Writing the thermal velocity 
vt in terms of Mach number M we obtain: 



1.44M 3/4 . 



(52) 



This formula is in good agreement with the simulated com- 
pression ratio . It was found that it is valid also for smaller 
free paths, when the maximum energy of the high energy 
tail is larger. 

This compression ratio is close to the maximum possi- 
ble compression ratio of the idealized system that includes 
the infinitely thin gas sub-shock and the viscous precur- 
sor produced by energetic particles. Using the Ranke- 
Hugoniot conditions at the thermal sub-shock one can ob- 
tain for this case: 



r = 2.5 (M 2 /5) 



3/8 



1.37M 3/4 . 



(53) 



This number is also close to the value found in numerical 
simulations of shock waves, modified by the cosmic ray 
pressure (see e.g. Berezhko & Ellison fl999|) . The compres- 
sion ratio of the thermal sub-shock is 2.5. 

The comparison of our results and results obtained for 
this simplified description are shown in Fig. 6. The Eq. (14) 
was solved in the upstream region together with the Euler 
equations. The top panel shows the comparison of down- 
stream spectra obtained in these approaches. The velocity 
profiles are compared in the bottom panel. We find a very 
good agreement of these two methods. We estimate the 
injection efficiency at injection velocity Ui n j = 2u su b to 
be about 0.024. Here u su b is the sub-shock velocity. This 
number is in reasonable agreement with results of hybrid 
simulations (cf. Giacalone et al. 119931 Ellison et al. 1993) 
and the Earth bow shock observations (cf. Ellison et al. 
11550) . 

We have also performed modeling of a quasi-parallel 
shock. The results obtained for a Mach number of 7.75 
and an angle between the shock normal and magnetic field 
of 9 = 15° are shown in Fig. 7. 
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Fig. 6. Comparison of our model with the cosmic ray mod- 
ified shock approach for a parallel shock with Mach num- 
ber 7.75, T e /Ti = and energy independent free path 
A = 0.2. The downstream distributions for our model and 
for the cosmic ray approach are shown in the top panel 
by solid and dashed lines respectively. Plasma velocity u 
from our simulations, pressure of ions P and velocity pro- 
file for the cosmic ray approach are shown in the bottom 
panel by dotted, dashed and solid lines respectively. 
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Fig. 7. Simulation of the quasi-parallel shock with 9 = 
15°, Mach number M=7.75, T e /Ti — and energy inde- 
pendent free path A = 0.2. The forward and backward 
particle distributions N + and N- are shown in the top 
panel by solid and dashed lines respectively. Plasma ve- 
locity u, pressure of ions P and electric force F = qE\\ 
shown in the bottom panel by solid , dashed and dotted 
lines respectively. The total compression ratio r = 7.21 is 
obtained. 



For such oblique shocks the dip appears between the 
main part of the thermal distribution and the high-energy 
tail (see the top panel). The injection efficiency is the same 
as in the previous case. 15° is the maximal value of 9 that 
allows us to model steady state shocks with M = 7.75. 
Even this value is rather large, because it corresponds to 
the normal angle 40° upstream of the thermal subshock. 
For larger 9 this dip becomes more pronounced and the 
thermal subshock width drops to a grid step. In some cases 
some instability in the downstream region was observed. 
It is possible to model the larger values of 9 for the smaller 
Mach numbers (e.g. 9 = 30° for M = 3.87). Thus our ap- 
proach can be used only for quasiparallel shocks with the 



normal angle less than 40° just upstream of the thermal 
subshock. 

If the magnetic field is oblique enough, it prevents par- 
ticles from returning from downstream along magnetic 
lines and the heating of particles is impossible. Since it 
is necessary to have the energetic particles downstream 
for the existence of a shock, we expect that the transi- 
tion region will be squeezed to widths of the order of the 
ion gyroradius. Heating then may take place in a different 
regime, where the drift motion of particles perpendicular 
to magnetic field is essential. The injection efficiency also 
may be very different. 
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6. Conclusion 

Many astrophysical problems deal with collisionless 
plasma. Different low frequency instabilities may produce 
magnetic inhomogeneities that can scatter particles and 
play the role of thermal collisions. When unstable waves 
propagate preferably along the mean magnetic field, the 
resonant scattering of particles is weak in the vicinity of 90 
degrees pitch angle. In this case it is possible that particle 
distributions are almost isotropic in both hemispheres. 

Following the approach of Isenberg (|1997p we derive 
the general transport equation (12) (or equivalent Eq.(9)) 
which takes into account the motion of the medium and 
energy changes of particles. Our consideration is not lim- 
ited by high particle velocities and by small anisotropics 
of particle distribution. 

In the case of high particle velocities and small 
anisotropies the equation derived can be reduced to the 
standard cosmic ray transport equation. In the case of neg- 
ligible energy changes and advective transport the equa- 
tion may be transformed to the so-called telegraph equa- 
tion. 

The equations derived can be solved together with 
magnetohydrodynamic equations (see Sect. 3). We used 
the fluid approximation for the electron transport. This is 
justified if electrons are effectively scattered by magnetic 
inhomogeneities. If this is not the case, one can use the 
equations (12) for electron transport. In this way electron 
injection at collisionless shocks can be investigated. 

The transport equations (12) with Eq. (19), (22) and 
(23) can be used for the solution of different astrophysi- 
cal problems when the self-consistent determination of the 
flow velocity, magnetic field etc. is necessary. 

We apply the equation derived to investigate particle 
acceleration and injections at astrophysical shocks. For 
the prescribed velocity profile we have found the exact 
analytical solution (see Sect. 4). 

We also have performed the modeling of collisionless 
shocks and found a formula for the shock compression 
ratio (see Eq.(52)). This ratio is very close to the maximal 
possible value in the framework of the simplified approach 
used for cosmic ray modified shocks. The corresponding 
thermal sub-shock ratio is 2.5. 

We found that this feature is universal and does not de- 
pend on the maximum energy of the accelerated particles. 
This was checked in our simulations up to the maximum 
momentum p max ~ 100 mu\ which was limited by the step 
of our uniform spatial grid. This property will not change 
for higher maximum energies corresponding to the super- 
nova environment. On the other hand we found that the 
injection efficiency rj is not universal. For p max ~ 10 mui 
it is about 0.024 for an injection velocity equal to two 
sub-shock velocities. It slowly decreased as the maximum 
energy increases and is about 0.003 for young supernova 
remnants. This is because at small maximum energies the 
thermal ions "feel" a larger compression ratio than the 
thermal subshock compression ratio since the subshock 
and precursor width are not strongly different. 



We found that the equations derived can be used to in- 
vestigate injection and acceleration in quasiparallel shocks 
with the subshock normal angle 9 < 40°. The dissipa- 
tion at more oblique shocks is described by a different 
mechanism which takes into account the motion of parti- 
cles in the direction perpendicular to the magnetic field. 
The injection efficiency at such shocks is unknown, but 
it is possible that it is rather low. We confirm the pre- 
vious findings that quasiparallel shocks are very efficient 
accelerators with high injection efficiency (see e.g Ellison 
& Eichler 115531 Giacalone et al. 115551 Malkov & Drury 
2001) which does not depend on the angle 9 (Giacalone et 
al. 119971) . 
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Appendix A: Formal inclusion of Alfven heating 

The heating of the cosmic ray precursor is described by the 
following equations (Volk & McKenzie 1982, McKenzie & 
Volk HSU}: 



d n 
—up = 

oz 

du 
PU d~z 

dp, 



dP, 
dz 



dP c 
dz 



du 



dP r 



'dz 



dz 



(A.l) 



(A.2) 



(A.3) 



Here V a is the component of Alfven velocity parallel to 
the shock normal, P g and P c are the pressure of the gas 
and cosmic rays respectively, and 7 9 is the gas adiabatic 
index. The third equation describes the gas heating due 
to the damping of Alfven waves generated by the cosmic 
ray streaming instability. The second equation is the Euler 
equation of motion. 

The gradient of the cosmic ray pressure can be found 
from Eq. (A.2) and substituted into Eq. (A.3). For high 
Mach number shocks the solution can be written as 



Pn 



Pgl + PlUlVal 



la 



1 



la 



la 



la ~ 

1/2 



(A-4) 



Here P g \ and V a i are the gas pressure and Alfven velocity 
in the medium, in which the shock propagates with speed 

Ul. 

For high Mach number shocks the Alfven heating is 
essential only at the very beginning of the precursor. In 
the rest, the gas is heated adiabatically and the last term 
in Eq. (A. 4) can be neglected. The gas pressure P g \ and 
Alfven velocity V a \ can be expressed in terms of Alfven 
and sonic Mach numbers M a and M s respectively. As a 
result, we obtain the sonic Mach number M of the shock 
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without Alfvcn heating, which is equivalent to the shock 
with the Alfven heating: 



M 



M 2 s M a 



1 9 i 2 



(A.5) 



For 7 g = 5/3 we obtain the formula (50) from the main 
part. 

Appendix B: Numeric method 

We used the implicit finite difference scheme for the solu- 
tion of the one-dimensional version of Eq. (12): 



dn 



= L z n + L v n. 



(B.l) 



Here n describes the pair of forward and backward particle 
distributions n± = v' 3 N±, and L z and L v are the opera- 
tors in the coordinate and velocity space respectively: 



L z n = (u± ^v'b z | // ! . 



(B.2) 



, d ( F\\ ldu\ 
L v n = - V — \±- r ---jn ± + , T n T - v ± n±. (B.3) 

The finite difference version of Eq. (B.l) was solved in two 
steps (GodunovUnZIl): 



fc+1 / 2 k — rdif i,j , rd,if„i,j 
~ — n k+l/2 + - L z n k ' 



n k+l ~ n fe+l/2 



(B.4) 
(B.5) 



The parallel force Fi\ was recalculated at each time step 
according to Eq. (21). The medium velocity u\ at each 
time step was calculated according to the finite difference 
version of the Euler equation (19): 



Pk+l u k+l 



pl< _ pUM+i) 2 - pi+M'+i) 2 



Az 



pi — l 
r fc+1 



pi+1 



2Ai 



(B.6 



Here P\ is the sum of the ion and electron pressures. The 
ion pressure and plasma density were calculated from the 
particle distribution at the each time step. The electron 
pressure was found from the adiabatic equation of state. 
The Eq. (B.8) is the quadratic equation for the velocity 
u\ +1 . The gas pressure at the right boundary P\ was pre- 
scribed in accordance with the shock motion: 



pb _ pb 
r k — r k-l 



1 



(z s k - z, 



k-1) 



(B.9) 



Here is the position of the shock. When the shock 
moves to the right boundary, the pressure increases. This 
permits us to reach steady state. 

We use the hyperbolic tangent for the initial velocity 
profile. The initial particle distribution was Maxwellian 
with a temperature and density dependence correspond- 
ing to the momentum and mass conservation. At times 
of about 60-100 in dimensionless units the system reaches 
the steady state. We use 200 grid points in the z-direction, 
100 grid points for the logarithmic v' grid and the time 
step between 0.005 and 0.1 for different runs. We obtained 
the total momentum, energy and mass conservation with 
several percent accuracy. 



Here is the value of distribution n(z, v' , t) at z = iAz, References 



v = Vj, t = fcr, where r is the time step, Az is the 
grid size in z-direction and Vj is the value of velocity v' 
at the knot of the velocity grid with number j. We use 
the following finite difference analog of the operator 
L z n = -d/dzAn (Kota et al. 11552]) : 



6 k k 



2.J 



-A 



k a k - 



l\] : "l ' •') /A*, A\ < 0, 



(B.6) 



2 k k 



-<^)/A:. X Ml. 



(B.7) 



A similar operator was used for the finite difference analog 
Ly % * of the velocity operator (B.3). 

The numeric scheme (B.4), (B.5) with these operators 
has a third order accuracy on z and v' and first order 
accuracy on t. 
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